;; WRITE INITIAL CONDITIONS FOR MERCURY

openw, 1, "../init.txt", WIDTH=250

NT = 16384

printf, 1, [NT, 0.01, 0.1, 5.0, 0.0001]
printf, 1, [1,0,0,0,0,0,0]

N = 2048    ;; number of stars
Rin = 0.5  ;; inner edge of the disc
Rout = 1.5 ;; outer edge of the disc
hr = 0.01  ;; aspect ratio

m = 8.33333333333333d-06 ;; mass of a single star in BH masses
Rs = 3d-5                ;; "size" of stars

G = 1
R = ( (Rout-Rin)*randomu(seed,2*N) + Rin )
phi = 2.*!pi * randomu(seed,2*N)

x = R*cos(phi)
y = R*sin(phi)
z = R * hr * randomn(seed,2*N,/normal)


vk = sqrt(G/R)  ; mbh = 1

vx = -vk * (sin(phi) + hr * randomn(seed,2.*N,/normal)/3.)
vy =  vk * (cos(phi) + hr * randomn(seed,2.*N,/normal)/3.)
vz =  vk * hr * randomn(seed,2.*N,/normal)/3.

for i = 0, N-1 do begin
  printf, 1, m, x[i],  y[i],  z[i], vx[i], vy[i], vz[i]
endfor

Nm = NT - N+1
Mm = 0.000833

R1 = 1.5
R2 =  20

; first choose radius, following n propto r^-7/4, the probability
; distribution for r is then f(r) = 5/4 r^1/4 / (R2^5/4 - R1^5/4)

Rm = dblarr(Nm)  ; here will store the radii
i = 0

while Rm[Nm-1] eq 0 do begin
  rtry = R1 + (R2-R1)*randomu(seed,1)
  ytry = 1.25 * R2^0.25 / (R2^1.25 - R1^1.25) * randomu(seed,1)
  if ytry lt 1.25 * rtry^0.25 / (R2^1.25 - R1^1.25) then begin
     Rm[i] = rtry
     i = i + 1
  endif
endwhile

;; now decide the position of the bhs
phim = 2.*!pi * randomu(seed,Nm)
thetam = acos(2.*randomu(seed, Nm)-1.)

xm = Rm * cos(phim) * sin(thetam)
ym = Rm * sin(phim) * sin(thetam)
zm = Rm * cos(thetam)

;; Now choose the energy for each star, which distribute f(E) propto
;; E^1/4, one has to do this for each radius.  E is actually defined
;; positive, E = Gm/r - v^2/2

E = dblarr(Nm)-1

for i=0, Nm-1 do begin
  E1 = 0.      ;; minimum E, otherwise it'd be unbound
  E2 = G/Rm[i] ;; maximum E, with v^2 = 0 (assumes Mbh = 1)
;; now choose E randomly
  while E[i] lt 0 do begin
     Etry = E1 + (E2-E1)*randomu(seed,1)
     ytry = 1.25 * E2^0.25 / (E2^1.25 - E1^1.25) * randomu(seed,1)
     if ytry lt 1.25 * Etry^0.25 / (E2^1.25 - E1^1.25) then E[i] = Etry
  endwhile
endfor

;; with E given, calculate |v|
  vm = sqrt(2.*(G/Rm - E))

;; and give it a random direction (isotropic)
  phivm = 2.*!pi * randomu(seed,Nm)
  thetavm = acos(2.*randomu(seed, Nm)-1.)

  vxm = vm * cos(phivm) * sin(thetavm)
  vym = vm * sin(phivm) * sin(thetavm)
  vzm = vm * cos(thetavm)

for i = 0, Nm-1 do begin
  printf, 1, Mm, xm[i],  ym[i],  zm[i], vxm[i], vym[i], vzm[i]
endfor

close, 1
end
